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The spatial and orientational distribution in a dilute active suspension of non-Brownian 
run-and-tumble spherical swimmers confined between two planar hard walls is calculated 
theoretically. Using a kinetic model based on coupled bulk/surface probability density func¬ 
tions, we demonstrate the existence of a concentration wall boundary layer with thickness 
scaling with the run length, the absence of polarization throughout the channel, and the 
presence of sharp discontinuities in the bulk orientation distribution in the neighborhood 
of orientations parallel to the wall in the near-wall region. Our model is also applied to 
calculate the swim pressure in the system, which approaches the previously proposed ideal- 
gas behavior in wide channels but is found to decrease in narrow channels as a result of 
confinement. Monte-Carlo simulations are also performed for validation and show excellent 
quantitative agreement with our theoretical predictions. 


I. INTRODUCTION 

The propensity of confined self-propelled particles to accumulate at boundaries is a trademark 
of active matter and has been reported in many experiments on bacterial suspensions [1-3] as 
well as simulations based on various models [4-6]. Several disparate mechanisms have been pro¬ 
posed in explanation, including wall hydrodynamic interactions [1] and scattering due to collisions 
with the walls [7], though recent theoretical efforts have shown that the mere interplay of self¬ 
propulsion, stochastic processes and confinement is sufficient to explain accumulation [8-10]. With 
few exceptions, however, these models have necessitated particle diffusion, which in reality is nearly 
negligible in bacterial suspensions where stochasticity in the dynamics takes instead the form of 
run-and-tumble random walks [11]. 

Understanding the distribution of active particles in confinement is especially critical for de¬ 
termining the mechanical force per unit area exerted by the suspension on the boundaries, or 
so-called ‘swim pressure’. This novel concept, which has received much scrutiny recently, describes 
the entropic force that must be applied on containing osmotic walls to keep self-propelled particles 
confined. Models based on the virial theorem [12-14] and on direct calculations of the wall me¬ 
chanical pressure [15] in infinite or semi-infinite collections of spherical swimmers have all arrived 
at a simple ideal-gas law II,: for the swim pressure in the limit of infinite dilution: 

n,: = n(D t = nS-, 
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where n is the mean number density, £ is the viscous drag coefficient of a particle and Dt = Vq /3A 
is the long-time translational diffusivity of an unconfined run-and-tumble swimmer expressed in 
terms of its speed Vo and mean tumbling rate A [11]. Equation (1) and its extension to finite 
concentrations have proven useful to explain motility-induced phase separation in suspensions of 
self-propelled colloids [12, 16], though its general validity as a thermodynamic equation of state for 
the pressure of active matter remains controversial [17-19] and appears to be limited to unconfined 
spherical particles [13, 15, 20]. 

In this work, we analyze the simple case of a dilute suspension of athermal run-and-tumble 
spherical swimmers confined between two parallel flat plates. We propose in §11 a kinetic model 
based on two probability density functions describing the spatial and orientational distribution 
of the particles inside the gap and at the walls, which are coupled via flux conditions and only 
account for the effects of swimming and orientation decorrelation by tumbling. Further, our model 
implicitly captures hard-wall steric interactions without requiring the use of a soft potential to 
describe wall collisions as in previous theories [15, 20]. A semi-analytical solution method is outlined 
in §111, which provides the full probability density functions and allows for a direct calculation of 
the mechanical swim pressure exerted on the walls in terms of the polarization of the surface 
distributions. Results for the distributions and swim pressure are presented in §IV, where they are 
shown to compare very favorably with Monte-Carlo simulations. 


II. PROBLEM DEFINITION AND THEORETICAL MODEL 

A. Problem formulation 

As a minimal model for an active suspension in confinement, we consider a dilute collection 
of self-propelled spherical particles confined between two infinite parallel plates separated by a 
distance 2 H (see figure 1). The swimmers are non-Brownian and simply perform a run-and- 
tumble random walk: straight runs of duration r at constant velocity Vo along the unit director p 
alternate with instantaneous tumbling events causing random and uncorrelated reorientations of 
p. The time r between tumbles is an exponentially distributed random variate with mean A -1 , 
where the tumbling rate A is assumed to be independent of position and orientation. To elucidate 
the interplay between run-and-tumble dynamics and confinement, we focus on the dilute limit and 
entirely neglect interparticle interactions. Particle-wall interactions are purely steric: as a swimmer 
meets one of the two surfaces, the normal component of its swimming motion is cancelled by a 
hard-core repulsive force causing it to stay at and push against the wall until a subsequent tumbling 
event reorients it into the bulk. Tumbling events occurring at the walls can lead to reorientation 
into the wall or into the bulk, so that a particle at a surface may need to undergo several tumbles 
before it is able to escape. 

There are only two length scales in the problem: the mean run length £ r = Vo A -1 and the 
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FIG. 1. Problem definition: run-and-tumble particles are confined between two flat plates separated by 2 H. 
The distribution of particles is a function of z and q = p • z = cos 6 G (—1,1). Orientations pointing towards 
the top and bottom walls are parametrized by = q and q I = —q, respectively, both defined in (0,1). 

channel width 2 H. We define their ratio as the Peclet number Pe = l r /2H = Vq/2\H, where 
the two limits Pe — > 0 and Pe —> oo describe weak and strong confinement, respectively. Due 
to the symmetry of the problem, the distribution of particles in the channel only depends on two 
degrees of freedom: the wall-normal coordinate z G (— H, H) and the wall-normal component of 
the particle director q = p • z = cos 8 G (—1,1). It is convenient to distinguish particles pointing 
towards the top and bottom walls, and to this end we divide the unit sphere of orientations into two 
hemispheres and define two distinct orientation coordinates c/ = q G (0,1) and = —q G (0,1) 
on each hemisphere for particles pointing up or down, respectively, as depicted in figure 1. 

The distribution of particles in the channel is then fully described by a bulk probability density 
function z , q) and by two surface probability density functions 'ifjl(q^) and ipi(q^), which are only 
defined over half of the orientations since the surfaces cannot sustain a concentration of particles 
pointing towards the bulk. By symmetry, we expect 

i/>(z,-q)=il’(-z,q), 'ip{z,q t ) =ip{-z,q l ) and (2) 

for (/l = r/1. Next, we describe the coupled bulk/surface conservation equations satisfied by these 
distributions, together with the appropriate boundary conditions. 

B. Bulk conservation equation 

The steady bulk probability density function ^( 2 :, q) satisfies the conservation equation 

V o <1 Ob <?) = -A tp (z, q) + ^ J ^ A ip {z, q) d q. (3) 

The left-hand side describes transport along z due to self-propulsion. Run-and-tumble dynamics 
is captured by the right-hand side, where the first term accounts for depletion due to swimmers 
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tumbling away from orientation q. and the second term for restoration due to swimmers tumbling 
from orientations q' into q. It is also useful to define the orientational moments of order j of 
the bulk probability density function on the full sphere and on the upper/lower hemispheres of 
orientations as 

Mj(z) = j qi i/i (z, q) dq and = J (q^y ip(z, q^) d^, (4) 

and we note that the zeroth, first and second moments correspond to the concentration, polariza¬ 
tion, and nematic order parameter fields: 


c(z) = M 0 (z), 

c n (z) = Ml\z), 


m(z) = M\(z), 
mJ^(z) = M\^{z), 


S(z) = M 2 {z), 
S n (z) = Ml\z). 


(5) 

( 6 ) 


By symmetry, it is straightforward to see that full moments of even order are even functions of z 
whereas those of odd order are odd functions. With these notations, the bulk conservation equation 
(3) simplifies to 



i!)(z,q) 


-ip{z,q) + ±c(z). 


(7) 


C. Surface conservation equations 


Similarly, conservation equations for the steady surface probability density functions at the 
walls can be written. We first define the surface concentration and polarization as 


c s = f and m s = f q^ (q^) dq^, 

Jo Jo 


( 8 ) 


and note that the values of c s and rn s are the same at both walls. With these notations, the 
conservation equation at the upper wall (z = +H) reads 


y 0 <? t = a - 2 c s 


(9) 


and a similar equation holds at z = —H. The right-hand side in equation (9) describes tumbling 
processes at the wall. The left-hand side, on the other hand, captures the flux of particles that enter 
the surface from the bulk by self-propulsion, and is therefore proportional to the bulk probability 
density function ^{H,q^) next to the wall. Evaluating the zeroth and first orientational moments 
of equation (9) yields simple relations between c s and m s and the values of the bulk moments in 
the vicinity of the wall: 


c s = 2 £ r mJ(H), 


m s 


\rrJ{H) + S\H) 


( 10 ) 
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D. Boundary condition and particle number conservation 

Equation (9) can be interpreted as a boundary condition for orientations pointing into the 
wall. For orientations pointing away from the wall, the swimming flux away from the wall must 
be balanced by tumbling of particles from the surface towards the bulk. Simply stated, particles 
on the surface that tumble to an orientation pointing into the bulk are transported away by self¬ 
propulsion. This leads to the additional condition 

= \\c s or l r q l i’{H 1 q l ) = \c s . (11) 

As c s is constant and finite, this condition suggests divergence and discontinuity of the bulk prob¬ 
ability density function for orientations parallel to the wall (r/1 —y 0), as will indeed be verified in 
our analytical solution and stochastic simulations. 

Finally, the above system of equations for the bulk and surface distributions is supplemented 
by a constraint on the total number of particles in the channel: 

2 c s + [ c(z)dz = N, (12) 

J-H 

where N is the total particle number in a vertical slice of unit horizontal cross-section. 


III. METHOD OF SOLUTION AND SWIM PRESSURE CALCULATION 
A. Integral equation for the moments 

We now outline a solution method for the system described in §11. As a first step, we derive 
an integral equation relating the bulk orientational moments to the concentration field. The bulk 
concentration equation (7) can be viewed as a linear inhomogeneous ordinary differential equation 
for ij)(z,q) where q is a parameter. We solve it by the method of variation of constants, treating 
orientations q' and q^ separately. After applying the boundary conditions (9) and (11), we obtain 
a general expression for the bulk probability density function: 

, . t i. c s \ (H±z) 1 , f z c(z') I" (z — z') 1 . , 

>=2W* exp rTTJ-j d *- ■ < i3 > 

Note that the bulk and surface concentrations c(z) and c s both appear on the right-hand side 
and are still unknown. However, equation (13) shows that their knowledge entirely specifies the 
bulk distribution ip(z,q). The bulk moments of order j on both hemispheres of orientations are 
immediately obtained by integration: 


U„2wi exi r 


l r q U \ 


= ^ + i 


H±z 


± 1^1 a/, 


where £j is the exponential integral function defined as 


£j(z) = 


‘ exp - J d u. 
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Finally, the moment of order j on the full sphere of orientations can be shown to be 


Mj( z ) — ( £j +1 


H + z 


+ 8 


7+1 


H-z 


+ 


rH 


'—H 


c(z') 
2 L 
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7+1 


z — z 


d z'. 


(16) 


B. Bulk concentration profile 


Setting j = 0 in equation (16) immediately provides an integral equation for the yet unknown 
concentration profile: 


/ \ C-s 

c(z) = w r 


Si 


H + z 

L 


+ Si 


H-z 


+ 


r H 


l-H 


c(z 0 
2 L 


Si 


z — z 


d z'. 


(17) 


Dividing through by c s , we obtain an equation for c(z)/c s that can be solved numerically. For finite 
£ r , we find that an approximate solution is easily obtained iteratively by casting equation (17) in 
the form Ck+i(z)/c s = f[ck(z)/c s \, starting with an initial guess which we take to be cq{z) = 0. In 
strong confinement (large Pe), the solution converges in 0(20) iterations, though more iterations 
are required in wider channels. 


C. Surface concentration 

To complete the solution, the value of the surface concentration c s must be calculated. To this 
end, we make use of a crucial property of the system, namely the overall isotropy of the suspension. 
Indeed, the spatially averaged orientation distribution Q(q) must be isotropic as reorientation due 
to tumbling is completely uncorrelated and is unaffected by the presence of the walls. This is 
expressed mathematically as 

/ H 

^ip(z,q)dz = (18) 

which can be combined with the surface conservation equation (9) to provide an equation for c s . 
The solution to the problem then proceeds as follows. Solving equation (17) using the iterative 
procedure outlined above provides a solution for c(z)/c s . This can be inserted in equation (13) 
to obtain ip(z, q) /c s , which can then be substituted into the overall isotropy condition (18) to 
solve for c s . As a final step, the surface probability density function ip s can be determined using 
equation (9). Solutions obtained by this method are presented in §IV, where excellent agreement 
with results from Monte-Carlo simulations will be shown. 

D. Swim pressure calculation 

The above formulation provides a direct way of estimating the swim pressure in the system, 
which is simply the force per unit area exerted by the particles at the walls as they push on the 
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surface. Specifically, the normal component of the motion of each particle at the upper wall is 
resisted by a force where £ is the viscous drag coefficient of one particle [12]. Knowing the 

surface probability density function ijj ' s , an expression for the swim pressure is then easily found as 

IK= f 1 CV 0 qW 8 tf)dql =(V 0 m s , (19) 

Jo 

where m s is the surface polarization. Using equation (10), this is also expressed in terms of bulk 
variables as 

n s = CU 0 lr \\m n {±H) + S n (±i7)l = C*f \\m n (±H) + 5 n (±i7)j . (20) 

J A L 

In bulk unconfined systems, previous models have led to the ideal-gas pressure 11,; of equation (1), 
which contains no information on particle orientations due to isotropy but follows the same scaling 
as equation (20). To compare both predictions, we define a dimensionless pressure as the ratio of 
equations (20) and (1): 

p = 5 = ^ = ^[^ (±ff)+5ti(±ff) ]' <21) 

where n = N/2H is the mean number density in our system. V — 1 quantifies the departure from 
the ideal-gas swim pressure. We will see in §4 that V —>• 1 in very wide channels {Pe —> 0), but 
deviates from 1 when Pe > 0 as a result of confinement. 

IV. RESULTS AND COMPARISON TO SIMULATIONS 

A. Simulation method 

To validate our model, we also perform Markov-chain Monte-Carlo simulations of run-and- 
tumble swimmers between two hard walls. During a run of duration r, the swimmer trajectory 
simply evolves as x{t+At) = x{t) + Vop At where At is a short time step. Each run is then followed 
by a tumbling event, where the new orientation vector p is picked randomly on the unit sphere. The 
time r between two consecutive tumbles is drawn from an exponential distribution with cumulative 
distribution function F{t) = 1 — exp [—At]. When a swimmer meets a wall, it remains there and 
continues to tumble until it reorients towards the bulk and swims away. Time-averaged bulk and 
surface probability density functions were extracted from orientational and spatial histograms, and 
convergence was checked with respect to At and to the duration of the simulation. 

B. Theoretical and numerical results 

Solutions for the bulk concentration profile are depicted in figure 2, where both the full con¬ 
centration c(z) and the partial ‘up’ concentration c^(z) are plotted for various values of the Peclet 
number, which measures the degree of confinement. The full concentration profiles in figure 2(a) 
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z/H z/H 


FIG. 2. Concentration profiles across the channel for various values of Pe = t r /2H: (a) full concentration 
c(z), and (b) partial ‘up’ concentration P(z). Solid lines show the semi-analytical solution of §111, and 
symbols are Monte-Carlo simulation results. 

show significant accumulation at the walls, with wall boundary layers whose thickness scales with 
i r . An interesting and unique feature of run-and-tumble particles is that accumulation occurs in the 
absence of polarization, and m(z) is found to be strictly zero throughout the channel (not shown). 
A non-zero polarization would indeed lead to a net flux of particles in the wall-normal direction, 
which cannot happen in a confined athermal system, unlike in Brownian suspensions where this 
flux can be balanced by diffusion [10]. In fact, averaging equation (3) over q immediately leads 
to the condition that m(z) = 0. The profiles also show the presence of a singularity in c(z) at 
the walls, which is a direct consequence of the boundary condition (11) and is also obvious from 
the solution (17) where £i(0) diverges. Concentration singularities were also predicted by Elgeti 
& Gompper [9], though their model did not capture orientation distributions. As confinement 
becomes significant and Pe increases, the bulk concentration decreases throughout the channel to 
reach nearly zero at Pe = 200, indicating that strongly confined particles spend most of their time 
at the boundaries. Excellent quantitative agreement is obtained between theory and Monte-Carlo 
simulations, thereby strongly validating our kinetic model. 

Figure 2(b) also shows the partial ‘up' concentration obtained by only counting particles pointing 
towards the top wall. The asymmetry of the profiles and the singularity at the bottom wall indicates 
that on average there are more particles pointing away from the wall than towards it inside the 
wall accumulation layers. However, in order to satisfy no net polarization in the bulk, this implies 
that those particles pointing towards the wall are more strongly polarized than those pointing 
away. This point is confirmed in figure 3 (a-b), showing the orientation distributions in the bulk in 
the vicinity of the top wall for orientations pointing away from and towards the wall. Figure 3(a) 
confirms the divergence of the bulk probability density in the neighborhood of orientations parallel 
to the wall (q' y —> 0) as expected from boundary condition (11), which is also captured by the 
simulations. The presence of this discontinuity can be rationalized as follows: particles that leave 
the surface at an orientation > 0 swim nearly parallel to the surface and therefore remain there 
much longer than particles leaving in other orientations. The distribution of particles pointing 
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FIG. 3. Bulk probability density at the top wall for (a) orientations pointing away from the wall and ( b ) 
orientations pointing towards it. (c) Surface probability density at the top wall as a function of q ( . Solid 
lines show the semi-analytical solution of §111, and symbols are Monte-Carlo simulation results. 


towards the wall in figure 3(6) shows no such singularity, but exhibits a finite peak at a critical 
value of whose origin remains unclear. The orientation distribution ipsio^) of particles on the 
top wall is shown in figure 3(c) and shows a preferential alignment normal to the wall rather than 
parallel to it. However, this distribution becomes nearly isotropic under very strong confinement 
(Pe = 1000), for reasons that we elucidate below. 

Taking moments of '(/4 (c/ ) provides the surface concentration c s and surface polarization rn s , 
which are plotted versus Peclet number in figure 4(o-6). Both quantities increase with increasing 
confinement, but asymptote as Pe —> oo. The asymptote for c s is N/ 2, meaning that in very 
narrow channels the particles spend all their time at the boundaries; indeed, the time 2H/Vq it 
takes them to cross the gap is infinitesimal compared to the mean run time A -1 . This is also 
consistent with the decrease in the bulk concentration seen in figure 2(a). In this limit, particles 
tumbling away from one wall reach the other wall nearly instantaneously, leading to an isotropic 
surface orientation distribution in agreement with figure 3(c), hence the asymptote of iV/4 for the 
wall polarization. 

Lastly, the dependence of the dimensionless swim pressure V on the degree of confinement is 
illustrated in figure 4(c). In the limit of weak confinement (H 3> £ r or Pe —>• 0), the swim pressure 
is seen to tend to the ideal-gas law of equation (1) in both our model and simulations: V —>• 1 or 
n s — > n.j. This corresponds to the limit of a single wall where the gap width H plays no role, and 
validates the results of previous studies in infinite or semi-infinite systems for which the expression 
for n* was first derived [12, 15]. Confinement, however, causes a decrease in the swim pressure, 
which in fact tends to zero for fixed n in very narrow gaps. The high-Pe asymptote for m s describes 
the limiting behavior: 


V 




-1 


i.e. lie —} 


3 „ _ lTT nH(V 0 N(V 0 
4 Pe Ul = ^ = — 


( 22 ) 


as Pe —> oo (or H —> 0), which corresponds to IV/2 particles pushing with an average force of 
QVq/2 against each wall. The decrease in pressure and the details of the asymptote agree with the 
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FIG. 4. (a) Surface concentration c s , (b) surface polarization m s , and (c) dimensionless pressure V as 
functions of Peclet number Pe = t r /2H. Solid lines show the semi-analytical solution of §111, and symbols 
are Monte-Carlo simulation results. 


previous two-dimensional results of Yang et al. [13], who also verified them in numerical simulations 
of self-propelled disks. They are also consistent with the study of Ray et al. [17], who analyzed 
the force on two nearby parallel plates in an active particle bath and proposed that the pressure 
inside the gap in a one-dimensional system with constant run length goes as II,/(1 + Pe). 


C. Summary and discussion 

We have presented a simple continuum model for a dilute suspension of spherical run-and-tumble 
particles confined between two hard walls and interacting via purely steric forces with the walls. 
The model improves upon our previous theory for confined Brownian suspensions [10] by allowing 
us to address the limit of zero temperature for the first time within a continuum framework and 
by incorporating a more realistic treatment of surface interactions and exchange processes between 
surfaces and the bulk without the need for a soft potential [15]. This description also provides 
a direct and simple way of calculating the mechanical swim pressure exerted on the walls. We 
have outlined an elegant approach to derive a semi-analytical solution for the probability density 
functions, and demonstrated excellent quantitative agreement between our model and results from 
discrete Monte-Carlo simulations. 

Our theoretical predictions and simulation results have highlighted several striking features 
of confined suspensions of run-and-tumble particles, namely the presence of a singularity and 
discontinuity in the bulk probability density function for orientations nearly parallel to the walls 
in the near-wall region, and the existence of a concentration boundary layer of thickness of the 
order of i r that actually diverges at the walls. Our pressure calculations were shown to match 
the recently proposed ideal-gas equation of state of active matter in wide channels, thus further 
validating this ideal-gas law and confirming the prediction that the precise nature of particle- 
wall steric interactions has no impact on the wall mechanical pressure for spherical particles [15]. 
We demonstrated, however, that confinement leads to departures from this ideal behavior and 
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specifically to a decrease in the swim pressure, which in fact vanishes in the limit of an infinitely 
narrow gap. In this case, we found that swimmers spend all their time at the boundaries, which 
provides the basis for previous models of strongly confined systems that only account for the surface 
distribution of swimmers [21]. 

While capturing the salient features of confined active suspensions, the problem under consid¬ 
eration remained minimal. Yet, the kinetic model presented here could be further modified to 
incorporate other effects and provide a more realistic description of biological or synthetic active 
systems. In particular, many active particles are rod-shaped and therefore also incur an aligning 
torque as they interact with boundaries. Recent theoretical work has shown that the wall pressure 
is modified in that case and becomes dependent upon the precise nature of particle-wall interac¬ 
tions [20]. In addition, experiments show that the surface-to-bulk tumbling of biological swimmers 
as well as certain types of synthetic swimmers is not uncorrelated but rather results in the prefer¬ 
ential release of the particles near a specific angle [22, 23]. Incorporating such details in our model 
is straightforward and would modify the distribution of particles near the walls with unexpected 
consequences for the mechanical pressure. Our basic model, validated here in the dilute limit, 
could also be modified to account for hydrodynamic couplings and to study the structure of the 
self-generated flows and collective dynamics of interacting active particles in confinement. Extend¬ 
ing the model to non-planar boundaries, whether concave or convex, is not as straightforward but 
would be of great interest for the theoretical description of active particle transport in complex 
geometries or of their interaction with and transport of passive payloads. This rich avenue is the 
focus of our current work. 
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